Data Analysis for Emma’s Cardiac Function Data

By Chance Yan on 4/2/25

This document is to show how I processed the cardiac function data collected by Emma Rawson.

  1. Data is given to me in a raw form. Essentially it’s a tdms file that we can’t read.

  2. We have to convert this tdms file into a txt file that we can read into an R script. To do this I used a macro script on Microsoft excel (convert tdms to text.xlsm).

  3. After we have our txt file (which is just a large list with a bunch of number readings) we can use the heart rate analysis script that Chris had created. This is the step where we determine what counts as a heart rate or note. To make things non-biased, I was sure to not look at the meta data before analyzing the heart rates and I also looked for repeated patterns that represented HRs. Some of the data had a lot of noise to it (making the signal unclear) in which case I either used a noise suppression tool to block out some noise or threw the signal out.

  4. After dictating what was a HR or not, I then pushed the data through a temperature matching script. I do this because when I highlight heart rate from an individual, I’m doing so from 10 minute intervals. I don’t highlight the whole 10 minutes (because that would give a poor signal), but only a section. By using the temperature matching script I’m able to tell where the HR begins exactly and where it ends leading to a higher accuracy in the final product (thus matches the HR to the temperature).

  5. After collecting the HR data and temperatures that matched with them, I can begin using the data into a readable form and data analysis (this script). Below is the start of that. I begin with combining the data and making sure I have everything into a data frame.

library(ggplot2)
library(patchwork)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(dplyr)
library(ggrepel) 

#combining data
{
  setwd("/Users/chanceyan/Documents/R/EmmaHR/EmmaHR/New_CardiacF(x)_converted")
  getwd()
  
  #creates data frame with all fulldata (heartrates of all trials), doesn't include trial 30 currently
  fulldata <- list.files(pattern = "fulldata.csv") 
  heartrates <- do.call(rbind, lapply(fulldata, read.csv))
  
  #loads and assigns ID to each invidividual
  require("readxl")
  setwd("/Users/chanceyan/Documents/R/EmmaHR/EmmaHR")
  my_data <- read_excel("Combined.xlsx")

  #removes all discards
  heartrates=subset(heartrates, confidence !="Discard")
  
  for(i in 1:nrow(heartrates)){
    if(heartrates[i,9] == "Flatline"){
      heartrates[i,6] = 0
    }
  }
  
  #assigns ID to each individual in the heartrates data frame
  heartrates$snail_id = paste(heartrates$Input, heartrates$Trial, sep="_")
  
  #combine heartrates df and meta df
  heartrates=merge(heartrates, my_data, by="snail_id", all.x = TRUE) 
  
  #heartrates$ind = NA
  individuallist <-unique(heartrates$snail_id)
  heartratesbackup <- heartrates
}
## Loading required package: readxl

There are 96 total individuals that Emma tested. 24 in each group and treatment.

a <- subset(heartrates, heartrates$confidence == "Flatline" | heartrates$confidence == "Zero")
table(a$population, a$group)
##     
##      Chronic Reversal
##   NC      24       24
##   NH      24       24

After I have put everything into the heart rates data frame, I remove unwanted signals (signals that were poor quality or did not run well with the script (mentioned later)). Looking at the groups again, we see there are less individuals, but still a large amount per group and treatment. After getting rid of the poor quality ones we see there are 85 signals.

{
  heartrates <- heartratesbackup
  #removal of problem ones that wouldn't go through code
  heartrates = subset(heartrates, snail_id != "0_Trial12")
  heartrates = subset(heartrates, snail_id != "1_Trial8")
  
  #Runs, but bad data/shape (noted poor signal)
  heartrates = subset(heartrates, snail_id != "3_Trial3")
  heartrates = subset(heartrates, snail_id != "5_Trial2")
  heartrates = subset(heartrates, snail_id != "1_Trial6")
  heartrates = subset(heartrates, snail_id != "4_Trial12")
  heartrates = subset(heartrates, snail_id != "5_Trial12")
  heartrates = subset(heartrates, snail_id != "4_Trial16")

  #Decreasing from the start (perhaps double check or add more points?)
  heartrates = subset(heartrates, snail_id != "0_Trial4")
  heartrates = subset(heartrates, snail_id != "0_Trial14")
  heartrates = subset(heartrates, snail_id != "2_Trial4") #maybe add more points and it'll be okay
  
  individuallist <-unique(heartrates$snail_id)
}

a <- subset(heartrates, heartrates$confidence == "Flatline" | heartrates$confidence == "Zero")
table(a$population, a$group)
##     
##      Chronic Reversal
##   NC      21       23
##   NH      19       22

I will translate the data to ABT since it will make the results cleaner and since the script might only work for ABT values…

#translates to ABTs
library(weathermetrics) #Install this package if you dont have it already
heartrates$x=1000/celsius.to.kelvin(heartrates$temperature) #creates x axis with temp in 1000/kelvin
heartrates$y=log(heartrates$Heartrate+1) # creates y axis with nat log of hr
individuallist <-unique(heartrates$snail_id)

These 3 graphs are meant to graph out all individuals HRs. The first one is for ABT individuals and mapping out where the break points would be without script limitations (meaning they can have 1-4 break points, but up to the script discretion). The second one would be the same but for regular data. The third one would be for ABT, but here you are defining how many breaking points YOU think are in the signal. I think in Chris’ work he only had individuals with one or two breaking points (so only one).

#selgemented function (ABT)
library(segmented)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:patchwork':
## 
##     area
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
#for(i in individuallist[1:length(individuallist)]) {
#  indiv=subset(heartrates,heartrates$snail_id==i)
#  out.lm<-lm(Heartrate~temperature, data=indiv)
#  print("-----------------------------------------") #I added these "print" lines to create annotations in the console, but they can be removed
#  print(i)
#  print("
#        ")
#  os<-selgmented(out.lm, Kmax=3, type="bic", bonferroni=TRUE)
#  plot(Heartrate~temperature, data=indiv, main=i)
#  plot(os, add=T, lwd=1.6)
#  points(os, add=T, lwd=2,col=c("red","blue","green","yellow"),cex=3)
#}

#selgemented function (regulatr data)
#library(segmented)
#for(i in individuallist[1:length(individuallist)]) {
#  indiv=subset(heartrates,heartrates$snail_id==i)
#  out.lm<-lm(Heartrate~temperature, data=indiv)
#  print("-----------------------------------------") #I added these "print" lines to create annotations in the console, but they can be removed
#  print(i)
#  print("
#        ")
#  os<-selgmented(out.lm, Kmax=4, type="bic", bonferroni=TRUE)
#  plot(Heartrate~temperature, data=indiv, main=i)
#  plot(os, add=T, lwd=1.6)
#  points(os, add=T, lwd=2,col=c("red","blue","green","yellow"),cex=3)
#}

#In the below script, the "segmented' function fits a different number of breakpoints for each individual as specified in an "n.breaks" column on the datasheet
#use this once you have determined the number of breakpoints for each individual, and want to double check you picked the right amount
for(i in individuallist[1:length(individuallist)]) {
  indiv=subset(heartrates,heartrates$snail_id==i)
  out.lm<-lm(Heartrate~temperature, data=indiv)
  os<-segmented(out.lm, npsi=first(1))
  plot(Heartrate~temperature,data=indiv, main=first(indiv$snail_id))
  plot(os, add=T, lwd=1.6)
  points(os, add=T, lwd=2,col=c("red","blue","green","yellow"),cex=3)
}
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

The following code is responsible for obtaining the break point values received from the graphs above and putting them into the data frame. There are also other values you receive as well such as slope, but for the analysis right now it seems we’ll only use breaking point temperature.

#collect data from regression lines huge code from here
{
  library(weathermetrics)
  
  heartrates$snail_id=heartrates$snail_id
  
  #divide your data into three different files based on the number of breakpoints (specified in the column "n.breaks.arrhenius")
  singleonly=subset(heartrates,Breakpoints=="1")
  
  #Create seperate individual lists for each file 
  singleonly <-unique(singleonly$snail_id)
  
  #script for single break individuals
  #---------------
  
  # 3. run loop - subsets the data frame for each sample and adds the results to the lists for lm(), segmented(), 
  # predicted values from segmented lines, and breakpoints contained under $psi of the segmented results 
  
  #create/ clear storage files
  fit <- list()
  fit.seg <- list()
  predicted.seg <- list()
  breakpoints <- data.frame(Psi=as.numeric(),ABT=as.numeric(),SE=as.numeric())
  
  #script to calculate predicted heartrate at breakpoint, and slopes, in single bp individuals
  
  {breakpointhr <- list()
    endslope <- list()
    initialslope <- list()
    breakpointslist  <- list()
    intercepts <- list()
  } #creat lists to hold data
  
  for(i in singleonly[1:length(singleonly)]) {
    fit.seg.A1<-list()
    fit[[i]] <- lm(y~x, data=subset(heartrates,heartrates$snail_id==i))
    fit.seg[[i]]<-segmented(fit[[i]], npsi=1, data=heartrates, na.action=na.exclude)
    breakpoints[i,] <- fit.seg[[i]]$psi
    fit.seg.A1<-fit.seg[[i]]
    breakpointhr[[i]]= fit.seg.A1$coefficients[1] + fit.seg.A1$coefficients[2] * fit.seg.A1$psi[2] #gives y1)
    endslope[[i]]=fit.seg.A1$coefficients[2] # initial slope
    initialslope[[i]]=(fit.seg.A1$coefficients[2]+fit.seg.A1$coefficients[3]) # initial slope (need to change code for double endpoints)
  } #run script
  
  library (plyr)
  library(weathermetrics)
  {
    breakpointhr <- ldply (breakpointhr, data.frame)
    endslope <- ldply (endslope, data.frame)
    initialslope <- ldply (initialslope, data.frame)
    
    breakpoints$snail_id <- rownames(breakpoints) 
    breakpoints$snail_id=as.factor(breakpoints$snail_id)
    colnames(breakpointhr)[1] <- "snail_id"
    colnames(endslope)[1] <- "snail_id"
    colnames(endslope)[2] <- "endslope"
    colnames(initialslope)[1] <- "snail_id"
    colnames(initialslope)[2] <- "initialslope"
    colnames(breakpointhr)[2] <- "finalbphr"
    
    breakpointslist=merge(breakpoints, breakpointhr, by="snail_id", all = T)
    breakpointslist=merge(breakpointslist, endslope, by="snail_id", all = T)
    breakpointslist=merge(breakpointslist, initialslope, by="snail_id", all = T) # merge breakpoint data with breakpointhr and slopes
    breakpointslist$ABT=1000/breakpointslist$ABT
    breakpointslist$ABT=kelvin.to.celsius(breakpointslist$ABT) #converts predicted breakpoint back to celsius
    breakpointslist$finalbphr=exp(breakpointslist$finalbphr)
  }
  singlebreakpointlist=breakpointslist
  
}
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
#heartrates2 has all meta data backed up
#heartrates3 has just flatline rows of all individuals
heartrates = merge(heartrates, singlebreakpointlist,by="snail_id")
heartratesbackup <- heartrates

The next code chunk basically creates another data frame to obtain the flatline temperature and graphs initially looks at averages for ABTs and FLTs across groups and treatments.

#use the following code below when grapging Flatline graphs
flatlines = subset(heartrates, confidence == 'Flatline' | confidence == "Zero")

ggplot(data = heartrates, aes(x = population, y = ABT, color = Treatment)) + geom_boxplot() + facet_wrap(~group) + ggtitle("ABT differences across group and treatments")

ggplot(data = flatlines, aes(x = population, y = temperature, color = Treatment)) + geom_boxplot() + facet_wrap(~group)+ ggtitle("FLT differences across group and treatments")

table(flatlines$population, flatlines$group, flatlines$Treatment)
## , ,  = 22C
## 
##     
##      Chronic Reversal
##   NC      11       12
##   NH      12       12
## 
## , ,  = 27C
## 
##     
##      Chronic Reversal
##   NC      10       11
##   NH       7       10

Next Steps?

My note